Bottlenecks to vibrational energy flow in OCS: Structures and mechanisms 
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Finding the causes for the nonstatistical vibrational energy relaxation in the planar carbonyl 
sulfide (OCS) molecule is a longstanding problem in chemical physics: Not only is the relaxation 
incomplete long past the predicted statistical relaxation time, but it also consists of a sequence of 
abrupt transitions between long-lived regions of localized energy modes. We report on the phase 
space bottlenecks responsible for this slow and uneven vibrational energy flow in this Hamiltonian 
system with three degrees of freedom. They belong to a particular class of two-dimensional invariant 
tori which are organized around elliptic periodic orbits. We relate the trapping and transition 
mechanisms with the linear stability of these structures. 

PACS numbers: 34.30.+h, 34.10,+x, 82.20.Db, 82.20.Nk 



I. INTRODUCTION 

How docs vibrational energy travel in molecules? An- 
swering this question succinctly seems a hopeless task 
considering the complexity of interatomic interactions 
in a molecule. Yet even before scientists were bur- 
dened by this knowledge, the so-called statistical theo- 
ries posited the answer: Vibrational energy travels "very 
fast" and distributes itself statistically among the vibra- 
tional modes of a molecule, assumed to resemble an as- 
sembly of coupled oscillators, well before a reaction takes 
place. Reaction rate theories based on these assump- 
tions - known collectively as statistical or RRKM the- 
ories [l|, 0, 0, 3 ~ remain reliable working tools of the 
practicing chemist because they have been vindicated in 
an overwhelming number of chemical reactions. 

However, numerical studies of Hamiltonian systems 
have provided solid evidence 0, 0, 0, 0, 0, EH that the 
approach to equilibrium usually proceeds more slowly 
than predicted by statistical theories [ll|, E2| - and it 
is also nonuniform, showing intriguing fits and starts. In 
particular, for Hamiltonian systems with two degrees of 
freedom, the familiar picture of chaotic seas, rigid bound- 
aries in terms of noble tori [l3j . leaky barriers in terms 
of cantori [3, EH has been well-established in the liter- 
ature, and these structures are found to be the source of 
anomalous transport in such systems [H, E3] ■ 

Beyond two degrees of freedom, the transport picture 
in terms of phase space structures is less clear. However, 
the phase space of higher-dimensional systems shows sim- 
ilar features such as the abundance of periodic orbits, and 
a mixture of chaotic and regular regions, the latter being 
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characterized (under some hypothesis) by invariant tori 
of various dimensions. The KAM theorem [T^j states that 
these structures are in general robust with respect to an 
increase of the perturbation or equivalently to an increase 
of energy. Understanding transport properties has to rely 
on these robust structures which are encountered by any 
typical trajectory. Roughly speaking, the presence of so 
many periodic orbits explains why generic trajectories, 
even when the system is strongly chaotic, display long 
intervals of near-regular behavior alternating with fits of 
chaos-a hallmark of anomalous diffusion. 

The slow approach to equilibrium started to be ac- 
knowledged a little over fifty years ago with the investi- 
gation of the dynamics of coupled oscillators by Fermi, 
Pasta and Ulam who showed that the relaxationproblem 
is far more complex than anticipated 0, 0, 0, 0, 0, EH ■ 
In chemical physics, anomalous diffusion was first impli- 
cated in the intramolecular vibrational energy relaxation 
of the carbonyl sulfide OCS molecule [TTJ] . The numerical 
study of a classical Hamiltonian model of OCS shows very 
slow energy redistribution among the vibrational modes, 
even in the fully chaotic regime □J, disagreeing strongly 
with the fast timescales derived from traditional statisti- 
cal theory. The understanding of the dynamics was suc- 
cessfully achieved for a collinear model of OCS which has 
two degrees of freedom [SJSJM M, HI . However, se- 



vere technical difficulties 23,24 



251 ] have prevented such 



a level of understanding beyond two degrees of freedom, 
and in particular, for the planar OCS model, in which 
the molecule is allowed to bend. 

In this paper, wc analyze the dynamics of a model for 
the planar OCS which is a Hamiltonian system with three 
strongly coupled degrees of freedom. The aim is to iden- 
tify the relevant structures in phase space which are re- 
sponsible for trappings and escapes, strongly influencing 
the transport properties (most prominently, the redistri- 
bution of intramolecular energy among the three modes). 
For example, rapid diffusion through phase space takes 
place through the so-called accelerator modes [26|. In 
contrast, sticky structures [27| like resonant islands or 
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tori influence the dynamics by strongly slowing down 
the trajectories passing nearby. All these structures are 
responsible for anomalous diffusion and fractal kinetics 
in the system (for recent surveys, see Refs. [ljl [I?} and 
references therein). Identifying these structures and the 
mechanisms behind trapping, escape and roaming is es- 
sential for understanding the transport properties of a 
given system. Given that there are many such structures 
in a realistic system, the only realistic hope for forming a 
generally valid picture of transport is to locate invariant 
structures which are responsible for the main changes in 
the transport properties. 

The specific question we address is: What are the 
structures in the phase space of OCS acting as dynamical 
bottlenecks to the diffusion of chaotic trajectories? What 
are the structures allowing transitions to other parts of 
phase space? For three degree of freedom systems, these 
invariant structures can be invariant tori with dimen- 
sions zer o ( s tag nation points), one (periodic orbits), two 
or three [H, [23, Hil ■ They can also include the stable and 
unstable manifolds of these objects [3l|. How are invari- 
ant structures relevant in the phenomena of capture in 
chaotic systems? For planar OCS, we find that the bot- 
tlenecks and the transition mechanisms from trapped to 
hyperbolic behavior are provided by a particular class of 
two-dimensional tori and their unstable manifolds. These 
results were recently announced in a Letter [32j ■ 

The paper is organized as follows: In Sec. [Ill we briefly 
recall some basics of the Hamiltonian model for the pla- 
nar rotationless OCS molecule. We also summarize the 
main results obtained on the dynamics of OCS relevant to 
the transport properties (both in the planar and collinear 
cases). In Sec. IIII A| we illustrate the transitions which 
occur in the neighborhood of periodic orbits using sev- 
eral representations: Time series, time-frequency analy- 
sis, and Poincarc sections. The striking common feature 
exhibited by many trajectories support the idea of some 
kind of universal transition mechanism. In Sec. IIII Bl 
after summarizing our methodology, we investigate the 
neighboring phase space structures which strongly influ- 
ence the dynamics of these trajectories. 



II. THE OCS MODEL 

A. The Hamiltonian 

The dynamics of the planar model of carbonyl sul- 
fide (OCS) can be described by a Hamiltonian model 
with three degrees of freedom with three strongly cou- 
pled, non-separable modes: There are two stretching 
modes and one bending mode. Each mode is represented 
by a coordinate and momentum pair, which we define 
as: i?i = d(C, S) and Pi for the CS stretching mode, 
i?2 = d(C\ O) and P 2 for the CO stretching mode, and 
finally, a = Z(OCS) and P a for the bending mode of the 
molecule (see Fig.QJ. Hamiltonian model for the rota- 
tionless OCS molecule has been provided in [Tl|, H3- It 
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FIG. 1: Equipotential surfaces of the collinear configuration, 
given by V(Ri, R2, 71") = E [see Eq. ©]. From center out- 
wards, energies are E = 0.03, 0.06, 0.09, 0.10, 1.6, 1.75, 1.8, 
2.0. The energies studied in this article are 0.09 (below dis- 
sociation of the weakest bond) and 0.10 (above dissociation), 
and the corresponding equipotential contours are shown in 
bold. 



has the form 

H(R 1 ,R 2 ,a,Pi,P 2 ,P a )=T(R 1 ,R 2 ,a,P 1 ,P 2 ,P a ) 

+ V(R 1 ,R 2 ,a), ( ' 

where T is the kinetic energy and V is the potential en- 
ergy. The kinetic energy is quadratic in the momenta 
and is provided as 

T = '-jP? + Y p l + V3P1P2 cos a 

+ P 2 ( JH- ^ 2 ^jCQSoA 

a \ 2Rj + 2R\ R X R 2 ) 

p . (Pi , P 2 \ 
— M3 "a sin a I — + — I , 

where [ii are the reduced masses. Based on available 
experimental data, the analytic model for the potential 
energy surface has been proposed in [33j | . In summary, 
V is given by 

3 

V{R u R 2 ,a) =J2Vi(Ri) + Vi(R!,R 2 ,R 3 ) , (2) 

where Vi(Ri) are Morse potentials for each of the three 
interatomic distances R\, R 2 and R3 = d(S,0), and 

Vi{R) = A (1 - cxp [-pi(R - R°)]f . (3) 

Here, R® are the equilibrium interatomic distances, 
and i? 3 is given by R 3 (Ri, R 2 , a) = (R\ + R\ — 
2R\R 2 cos a) 1 / 2 . At equilibrium, the molecule is 
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collincar, therefore R® = i?° + Also, the interaction 
potential Vj assumes the Sorbie-Murrell form: 

3 

V I= A P(R U R 2 , R 3 ) l[ (1 tanh 7i [i^ — R°]) , 

i=l 

where P(Ri, i?2, i?3) is a quartic polynomial in each of 
its variables: 

P{R U R 2 ,R 3 ) = 1 + cf + eg 

+ C ijk R '. R j R k + c ^ll R i R j R kRl ■ 

All the coefficients (/ij, A, ft, 7i, A 4^' c g ' c ijL 
c ifLi) arc provided in Ref. We display the equipo- 

tcntial surfaces of V(R\, R2, Ri + R2) of the collincar 
configuration in Fig. [T] The equations of motion can be 
derived from Hamiltonian |T]) using the canonical Hamil- 
ton's equations. 



B. Summary of prior results on the OCS dynamics 

The classical models of both the collinear and the pla- 
nar (rotationless) carbonyl sulfide OCS molecule have 
been studied in detail in Refs. [IE 01, EE EE El EE El • 

The dynamics in the collinear configuration of OCS 
was first studied by Carter and Brumer [ill ]. They 
characterized the motion of this system at a number of 
energies, extending up to 20, 000 cm -1 (which amounts 
to E = 0.09 a.u.) A relaxation time, as defined in 
Refs. [13, EE EH, was estimated at 0.17 pico-seconds. 
However, after integrating trajectories for 2.4 picosec- 
onds, no relaxation to statistical equilibrium was ob- 
served. When this contradiction was investigated by 
integrating the equations for much longer times (up to 
45 picoseconds), two distinct timescales for relaxation 
were found, the longer of which characterized energy re- 
distribution that was incomplete even after 45 picosec- 
onds [35J. Even on the picosecond time scale, sudden 
transitions between relatively long-lived regions of local- 
ized mode energies were observed. Since this collinear 
model has two degrees of freedom, Davis and Wagner [35 1 
used Poincare surfaces of section as a visualizing tool 
for phase space structures. These revealed that even at 
high energy (E = 0.09), the system has a "divided phase 
space" , with coexisting regular and chaotic regions. They 
observed that trajectories can be trapped in restricted 
regions of phase space for many vibrational periods, af- 
ter which they would suddenly move to other regions of 
phase space to repeat the pattern. 

Progress came with the recognition that the then- 
recent lobe dynamics [II [H[ could help to explain 
non-statistical relaxation in two degree of freedom sys- 
tems [lH. When the strength of the perturbation (or 
equivalently, the total energy) is increased, the two- 
dimensional invariant tori of a Hamiltonian system with 
two degrees of freedom develop sets of "holes" with the 



systematics of Cantor sets. These holes, dubbed "can- 
tori" [3], form leaky barriers which can act as bottle- 
necks to phase space transport. These bottlenecks are as- 
sociated with broken tori with irrational frequency ratios, 
where those with "noble" number ratios being gcncrically 
the very last to be destroyed by an increasing perturba- 
tion (the supporting argument being that these numbers 
are the most poorly approximated by rationals [40lp. For 
OCS, their existence has been confirmed in Ref. [l8[ in 
a region between two resonances wco/^cs = 3/1 and 
wco/^cs = 5/2. The noblest irrational number between 
the rationals 5/2 and 3/1 is 2+7, where 7 = (s/E— 1)/2 is 
the golden mean [3, El and can be expressed as a con- 
tinued fraction of an infinite sequence of ones, also writ- 
ten as [1, 1, 1, 1, . . .]. These results obtained from classi- 
cal mechanical were confirmed using quantal wave packet 
calculation [4l|. However, these successful results could 
not be extended to the planar OCS due to severe tech- 
nical and computational difficulties [IE EiJ. Yet there 
were indications that this problem of intramolecular en- 
ergy flow in higher dimensions is also related to the reso- 
nant and non-resonant structures [2lL 36] . In particular, 
the relevance of Arnold's web in the diffusion of trajec- 
tories was highlighted. Among their conclusions are that 
transport is most rapid along low order resonance zones; 
transport is slow (diffusive) along high order resonances; 
it was conjectured that pairwise noble frequency ratios 
play a role of inhibiting transport along resonance lines. 



III. TRAPPINGS AND TRANSITIONS IN THE 
PLANAR OCS: BOTTLENECKS AND 
TRANSITION MECHANISMS 

A. Observations 

The complexity of transport processes in the collincar 
OCS model, revealed in the early investigations, suggests 
that a look into phase space structures such as periodic 
orbits or invariant tori is needed for a better understand- 
ing of these processes. Even if the measure of such invari- 
ant structures embedded into a chaotic sea is typically 
zero, the "neighborhoods" of influence around them can 
have relatively large measure and their finite-time prop- 
erties, as characterized by Lyapunov exponents, provide 
a quantitative picture of transport. The rationale goes as 
follows : An ensemble of trajectories, described by a den- 
sity function, which is centered in a finite volume around 
a periodic orbit, will evolve in finite time following this 
periodic orbit, and spreading predominantly in the direc- 
tion of unstable manifolds, exponentially in time with a 
rate equal to the local Lyapunov exponent. An orbit in 
the "neighborhood" of a periodic orbit, temporarily as- 
sumes or "shadows" the properties of this periodic orbit 
as a general consequence of dynamical continuity [42l ]. 
This temporary influence of periodic orbits can also be 
viewed as instantaneous time-periodic forcing, exerted by 
a periodic orbit. It is expected that, in general, typical 
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trajectories are trapped for longer times in the neighbor- 
hoods of linearly stable orbits. In what follows, we draw 
a dynamical picture of transport in OCS based on the 
determination of invariant structures in phase space and 
their linear stability properties. 



1. Density of periodic orbits 

A generic feature of Hamiltonian dynamics is the abun- 
dance of periodic orbits in phase space. Figure [21 rep- 
resents the averaged density in the configuration space 
(Ri, R2) of periodic orbit points on the Poincare section 
£ (defined in Sec. IIII A 3[) for planar OCS. A closer in- 
spection of this figure shows that the most prominent 
regions of stability surround short periodic orbits with el- 
liptic linear stability. A typical trajectory passes through 
this maze of periodic orbits, being trapped for some time 
according to local stability properties. The aim of this 
manuscript is to understand how a typical trajectory can 
be trapped and released locally around a given periodic 
orbit. In what follows we analyze the transport proper- 
tics in the neighborhood of an elliptic periodic orbit, like 
for instance O a , as shown by a circle in Fig. [3] 



2. Time-frequency analysis and stroboscopic mapping 

To examine the temporal features of trajectories 
we use time- frequency analysis [43j]. In what fol- 
lows, x designates a point in phase space, i.e. x = 
R2, a, Pi, P%, P a ). A finite segment of a trajectory 
can be represented by a sequence of phase space coordi- 
nates, {x„} n= i r . x„ = x(£ n ), visited by a trajectory 
at times t n . For a stroboscopic map, we take snapshots 
with a fixed time increment, t n+ \ = t n + A. It is nat- 
ural to scale the time increment A by the period T p of 
the organizing periodic orbit. We select A = T p /4. The 
time series of selected orbits are displayed in the bottom 
panels of Figs. [4] and [6] 

We study the instantaneous frequencies using wavelet 
decomposition. As described in Refs. [H, EH EH, the 
time-frequency analysis is based on a continuous wavelet 
transform of an observable f(t) 



Wf(t, 



f(rW 



dr. 



(4) 



We choose the mother wavelet tp, in the Morlet-Grossman 
form: ip(t) = e^'e"*' '' ,2a2 ' / '(cr 2 ^) 1 / 4 , with adjustable pa- 
rameters 77 and a. The time-frequency representation is 
obtained via a relation between the scale s and the fre- 
quency £ = 77/s. We consider the normalized scalogram 

Pwf(t,Z = ri/s) = \Wf(t,s)\ 2 /s, 

which can be interpreted as the energy density in the 
time-frequency plane. The ridges of Pw can be inter- 
preted as instantaneous frequencies, or more rigorously, 
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FIG. 2: Averaged density of periodic orbit points on the 
Poincare surface of section, projected onto the (Ri, i?2)-plane, 
weighted by the "local escape rate" jp , the sum of pos- 
itive Lyapunov exponents, A' p ' > or, in terms of Lya- 
punov multipliers for the periodic orbit p, given by 7+ = 



a 



|Af'|> 



where A^ p ' is an eigenvalue of DT^ 



evaluated at the periodic points. Periodic orbits with the 
following number of returns to the Poincare sections are de- 
termined: 1(4), 2(9), 3(10), 5(24), 7(26), 8(101), 11(40), 
13(33), 17(21), 19(43), 23(41), 29(34), 31(28), 37(43) where 
the number of orbits is shown in parentheses. Energy is set 
at E = 0.09. Lighter areas are dominated by more regular or- 
bits, darker by unstable orbits. The circle indicates the region 
located near O a where (7?i,J?2) ~ (3.6,2.3) where trapping 
and roaming is analyzed in Fig. [3] 
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FIG. 3: The periodic orbit O a at E = 0.09: projections in the 
(Ri, i?2)-plane (left panel) and in the (Ri,Pi), (i?2,P2) and 
(a, P«) planes (right panels). The blue curve in the (i?i,7?2) 
projection is the boundary of the energetically accessible re- 
gion. Dots indicate the location of the intersection with the 
Poincare section E. The periodic orbit O a is of elliptic-elliptic 
linear stability type (see Tab. U for details). 



the set of frequencies for a given time interval. In this 
section, two typical trajectories (whose initial conditions 
are specified in Tab. [J) which are initially close to el- 
ementary organizing periodic orbits, arc represented in 
Figs. 2] and [5] where the signal f(t) is chosen to be P\{t) 
or P2(i). It should be noticed that other choices of ob- 
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FIG. 4: Lower panel: Time series Ps(t), of a trajectory with 
initial coordinate ai(O a ) (see Tab.H]) near the periodic orbit 
O a . The energy is E = 0.09. Time is scaled to T^° a \ the pe- 
riod of O a - The integration time is T max = 512Tp 0<1 ' « 34 ps. 
Upper panel: Ridges of the time-frequency decomposition of 
P2(t). The frequencies of Pa(t) are denoted £p 2 , and are rep- 
resented in units of (Tp "') -1 . The shaded band locates the 
transition region. 
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FIG. 5: Two stable periodic orbits, Ob (upper panel) and O c 
(lower panel) for E = 0.1: projections in the (Ri , R2) -plane 
(left panel) and in the (Pi, Pi), (ife,^) and (a,P a ) planes 
(right panels) . These orbits are relevant in the trapping of tra- 
jectories with initial conditions 02{Ob) (see Tab.|TJ. The zero 
velocity curve (boundary in (Ri, R2) coordinates) is shown in 
blue. Dots indicate the location of the intersection with the 
Poincare section E. 



servables f(t) lead to the same qualitative features as the 
ones presented here, and in addition, these features are 
common to a wide set of other trajectories in the same 
neighborhood. 

Time-frequency analysis shows that each of these tra- 
jectories displays qualitatively distinct regions : some 
with approximately constant ridges in time, and others 
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"Note that in Ref. 1321 the second return map was considered so 
the stability indices of O a are half of the ones here. 

TABLE I: Initial conditions of the trajectories considered in 
the manuscript: The three periodic orbits O a , Ob and O c , 
and the two trajectories ai(O a ) and 02 (Oi,), one close to O a 
and the second one to Ob- First column: label of the initial 
conditions, value of energy E, and in case of periodic orbit, 
period T p and rotation numbers ui. Second column: Pi, R2 
and a. Third column: Pi, P2, and P a . 



with multiple and short ridges. These two regimes are 
clearly marked with transition intervals (highlighted by 
shaded bands in the figures). In Fig. 0] we observe a 
clear, sharp transition stage between trapped behavior 
(around O a ) and roaming behavior throughout a large 
portion of phase space. After some time spent around 
the periodic orbit, the trajectory seems to find an exit 
channel through a bottleneck. Gcncrically, any trajec- 
tory experiences multiple events of capture and escape 
(like the one in Fig. [5]). We have found that escape to 
the chaotic region proceeds in two stages, characterized 
by two different rates of escape. The transition inter- 
val is characterized by ttrap and t csc . The first ( "slow") 
stage, < t < ttrap, and the second ("fast") stage 
ttrap < t < t csc . The precise transition points located 
at ttrap and t esc may vary in different situations, but typ- 
ically t csc — ttrap <C ttrap ■ We notice that these trapping 
and transition stages although visible, were not as clearly 
apparent on the time series as on the time-frequency 
plots. 

In order to identify the phase space regions visited dur- 
ing the trapping and escape stages, we complement the 
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FIG. 6: Lower panel: Time series Pi(t) of a trajectory ini- 
tially close to Ob with initial coordinates <T2(Ob) (see Tab.|TJ. 
The energy is E = 0.1. The integration time is approximately 
512Tp° b \ Upper panel: Ridges of the time-frequency decom- 
position of Pi(t). We notice that they are three transition 
regions (located by shaded bands) in this case. 

time- frequency analysis by projections of segments of the 
time series in a two-dimensional plane given by two coor- 
dinates, e.g., (R U R 2 ), (Ri,Pi) or (R 2 ,P 2 ). InFig.0 two 
segments of the trajectories of Fig.S]are represented (left 
and middle panels), one corresponding to the trapped 
stage (to the left of the shaded band), and the other one 
to the trajectory after the escape process (to the right 
of the shaded band). It is shown that the trajectory is 
trapped into a small L-shaped region around the peri- 
odic orbit C a , and that after the transition point, the 
trajectory has access to a larger part of phase space with 
an apparent size of the order of the entire accessible re- 
gion. The same observation follows for the trajectory of 
Fig. [6l the stroboscopic plot of which appears in Fig. [8] 
By drawing tubes around it, we notice that the trajec- 
tory in both trapped segments sticks to particular regions 
around different periodic orbits. 



3. Poincare sections 

We use Poincare sections as another way to visualize 
multidimensional trajectories. Given that this Hamilto- 
nian system has three degrees of freedom, the Poincare 
section is four dimensional. We show below how two- 
dimensional projections of these sections can be used to 
gain insight into the dynamics (although this information 
is displayed less clearly for this system than for a system 
with two degrees of freedom). Given some scalar function 
U(x) of the phase space variables, we define this section 
£ to be the set of points x of a trajectory such that 

C/(x)=0, 

with x ■ dU/dx > 0. From two consecutive points 
x„ = x(i n ) and x„ + i = x(i„ + A(x„, t n )) on the Poincare 



section, we define a Poincare map Ty,, 

^s(x„) = x„+i. 

In what follows, we have used the surface £ defined by 

E/(x) = P a . (5) 

The argument for choosing this surface goes as follows: 
We study an energy range where the time series of the 
bending mode (a, P a ) oscillate around an instantaneous 
mean value (a) . Between each oscillation there is a turn- 
ing point where momentum P a vanishes. The only case 
which is not captured by the Poincare section is when the 
bending mode is "frozen" to a = tt and P a = 0, which 
corresponds to the collincar OCS. 

We choose a four-dimensional paramctrization of the 
surface of section E which consists of R\ , Pi, R 2 and P 2 ■ 
Setting P a = in Eq. ((6|), the equation 

H(R 1 ,R 2 ,a,P 1 ,P 2 ,0) = E, (6) 

with a constraint P a > is to be solved for 
a(R\,P\,R 2 ,P 2 ;E) numerically. There are two merits 
in using Poincare sections: First, representing projec- 
tions as a set of planar plots of canonically conjugated 
variables helps in perceiving the symplcctic symmetry of 
structures. Second, the section manifold is one dimension 
smaller than the energy manifold, and so are the maps 
of all the invariant structures. For instance, periodic 
orbits correspond to a finite set of points {x. n }n=i,...,N 
on S, and the dynamics visits these points in a cyclic 
manner, i.e. fs(x„) = x„ + i for n = 1, . . . , N — 1, with 
■^(xjv) = xi. Similarly two-dimensional tori correspond 
to closed curves on S. 

In Figs. [9l andflOl Poincare sections, projected on the 
planes (i?i,Pi) and (R 2 ,P 2 ), are drawn for the two tra- 
jectories considered in Figs. 2] and [5] These Poincare 
sections clearly show distinct one-dimensional curves 
(clearly visible in the insets of Figs. [9] and fTO]) in the 
transition stages (shaded bands on Figs. 0] and [5]). The 
tubes which we identify as two-dimensional invariant tori 
in phase space [Hj, represented on Figs. [7] and corre- 
spond to these one-dimensional curves (or more generally 
a set of one-dimensional curves) on the Poincare sections. 

In the trapping stages (around specific periodic orbits) 
like the ones in Figs. [7] and [8] (left panels), the rotation 
numbers are obtained from the frequency map analy- 
sis [ifl on the surface of section S. Dimensionless ra- 
tios of frequencies arise naturally in the Poincare map 
^-"ej and ratios of frequencies are called rotation numbers. 
The trapping stage can be characterized by a single ro- 
tation number (and its harmonics), implying that a two- 
dimensional torus is the relevant invariant structure in 
the trapping process. In the following, we determine such 
structures and highlight the family of two-dimensional 
tori which are relevant for the transport picture in this 
system. 
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FIG. 7: Left and middle panels: Stroboscopic plots of segments of the trajectory of Fig. [4] (for E = 0.09): Before the transition 
stage (left panel), and after (middle panel). Note how during the trapping (left panel) the chaotic orbit remains around the 
L-shaped stable periodic orbit O a (represented in red) before escaping into the chaotic zone (middle panel) . The trapping stage 
is inside the invariant structure (a two-dimensional torus) shown on the right panel. 




2.5 3 3.5 2.5 3 3.5 2.5 3 3.5 



2.5 












1 ~ 


R 2 

2 
















2.5 


3 


3.5 


2.5 3 3.5 

Ri 


2.5 


3 3.5 



FIG. 8: Left and middle panels: Stroboscopic plots of segments of the trajectory of Fig. [6] (for E — 0.1): trapped stage for 
t € [0, 110] (upper left), chaotic region for t £ [110, 350] (upper middle), trapped stage for t £ [350,460] (lower left), and chaotic 
region for t > 460 (lower middle). The trapping stages (around periodic orbits represented in red) are occurring inside the 
invariant structures (two-dimensional tori) shown on the right panels. 



B. Lower dimensional invariant tori 

1. A summary of the methodology 

It is well established that invariant structures in phase 
space play an important role in the transport proper- 
ties associated with Hamiltonian systems with two de- 
grees of freedom [l3|. In particular, the role of periodic 
orbits has been singled out in many experiments [47j |. 
Even if some aspects of this dynamical picture can be 
extended to systems with a larger number of degrees of 
freedom, it remains to address the role of invariant struc- 
tures which are not present in systems with one and two 
degrees of freedom, but are specific to three and more 



degrees of freedom. In three degree of freedom systems 
this new type of invariant structures takes the form of 
two-dimensional invariant tori. Observations described 
in Sec. IIII Al indicate that such tori close to elliptic peri- 
odic orbits play an important role. To have a qualitative 
description of dynamics near a periodic orbit, we con- 
sider a fixed point xo on the surface of section S, i.e. 
J'e(xo) = x o, corresponding to a point of a periodic or- 
bit. Near x , the Poincare map can be expanded into 
a linear part and a remainder: 

jF s (x) = JF s (x ) + J D^" s (x )(x - x ) + ft(x - x ) , (7) 

where DJ^fao) is the matrix of first order derivatives of 
the Poincare map, constrained to the surface of section 
and evaluated at Xq . All higher order terms in x — xo 
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FIG. 9: Poincare sections of the trajectory analyzed in Fig. [4] (for E = 0.09) on the (Ri, Pi)-plane (left panel) and on 
the (i?2, P2)-plane (right panel). It is apparent (in the insets) that before escaping to the external region, the trajectory is 
stuck in the neighborhood of a well-localized structure. The five points (in bold) are the intersections with E of a partially 
hyperbolic resonant periodic orbit which is responsible for the escape to the chaotic region through its unstable manifold. Two 
"tentacles" starting at the upper and lower parts of the structure are marked with broken lines connecting crosses to clarify 
what happens during the escape stage (shaded band in Fig. The resonant 2:5 periodic orbit has elliptic-hyperbolic stability 
with A = 0.113231998 per return (or A = 0.566159992 for the entire orbit) and a rotation number lu/ty = 0.35684077865194591 
per entire orbit. 
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FIG. 10: This figure illustrates the multiple capturing which is found generically for a randomly selected trajectory, as the one 
of Fig. [6] (for E = 0.1). The quasi-regular intervals are color coded: between iterations 1-47 (red) and 370-430 (blue). Note 
that the second interval draws two curves, because it takes two returns to the surface of section to draw this curve (i.e. it is an 
apparently connected curve for The two closed loops are the bottleneck torus of the bottom-right panel of Fig. [8] 



are collected in lZ(x — xo). Finite-time dynamics near the 
fixed point Xo are determined by the properties of the ma- 
trix DJ-'sfeo). Assuming that linearized approximation 
is effective, and discarding the remainder term from fur- 
ther discussions (the fully nonlinear problem with large 1Z 
is solved using the methodology outlined in Appendix IB"]) , 
we consider a closed curve f(s) on the Poincare section 
S defined on a torus sfT 1 , and consider the dynamics 
of x(s) = xo + (j(s) given by 

^s(x(s)) = x + eDFx(x )j(s). (8) 



If Ty, has at least one pair of eigenvalues in the form 
A = cxp(zttw), it is possible to find a 7(5) such that 
DJ-'ji'y(s) = j(s + lu € ) and \uj — ui t \ = o(e). Therefore the 
equation 

^■ E (x(«))=x(« + w) > (9) 

has a family of solutions, parametrized by the rota- 
tion number to. Equation ([9]) defines a torus as a loop 
on the surface of section £ with rotation number m. 
Even if D!F^(x ) has two pairs of eigenvalues of the 
form exp (±(,Wj) such an invariant loop close to Xq can 
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be found. More details on the determination of two- 
dimensional invariant tori are given in Appendix [Bl 



2. Invariant tori and their bifurcations: Bottlenecks 



In the cases discussed in Sec. IIII Al trajectories un- 
dergo a transition (after a trapping stage) in the vicinity 
of a nonresonant elliptic periodic point Xo, whether it 
is associated with O a or Ob- For each of these periodic 
points, the matrix Z3J r s(xo) has eigenvalues exp (±twi), 
exp(±tw 2 ) (numerical values are given in Tab. U) Pro- 
cesses associated with the escape from the trapping stage 
can be better understood by analyzing the tangent space 
of the elliptic periodic orbit O a that locally has the struc- 
ture of a direct product (center + center) T 1 x Ii x T 1 x I2 , 
with the periodic orbit at the origin. The elements of the 
two intervals Ii C 1 are rotation numbers w%, which are 
not unique in general: The choice is fixed by requiring 
lim^o^i = <<4\ where e is a measure of the "diameter" 
of the torus and uif are stability angles of the organizing 
periodic orbit. The Poincare map induces rotations 
on T 1 , r ui x 1 x r LU2 x 1, where r u is a rotation on T 1 
with the rotation number 10. Partial (or complete) res- 
onances are determined by one (or two) resonance con- 
ditions nui + muj2 + k = 0, where (n, m, k) are integers 
such that \n\ + \m\ + \k\ > 0. The most striking trap- 
ping effects are observed for partial resonances of the 
type T 1 x T x {0} x {0}, and {0} x {0} xT'x I 2 . They 
are two-dimensional manifolds (locally), and can be fo- 
liated by one-dimensional invariant closed curves, called 
hereafter "loops." We propose to investigate a resonance 
manifold by mapping out dynamical invariants that form 
its backbone structure. Choosing either of the two sit- 
uations, a resonance channel has been constructed by 
finding the two-dimensional invariant tori for u>i £ L . At 
a small distance from the periodic orbit we use infor- 
mation obtained from the linear normal form DJ--z(x.q). 
Once an initial loop is found, we follow the progress as 
the rotation number u> is varied continuously monitor- 
ing their stability properties. Local normal stability of 
each family of tori can be represented by plotting the 
maximal Lyapunov exponent A by solving the general- 
ized eigenvalue problem [See Appendix [Bl and Eq. (|B3|) ] . 
versus the rotation number uj. Such a plot for a family 
of two-dimensional tori, originating from C a , is shown 
in Fig. 1111 From Fig. 111! we obtain a transition point at 
ui/ir sa 0.481 in the form of a bifurcation of an invariant 
torus. Projections of two-dimensional invariant tori in 
the transition regions are shown in Figs. and [51 while 
corresponding loops in the surface of section E are shown 
in Figs. [5] and [TO] 

The transition stage (see Figs. [5] and [TUJ) indicates an 
exponential divergence resulting from an escape along the 
unstable branch of a hyperbolic manifold. The maximal 
Lyapunov exponent of the segment of a trajectory in the 
capture stage can be estimated by observing the dura- 
tion of capture, and the per-return Lyapunov exponent 
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FIG. 11: Maximal Lyapunov exponents of a family of two- 
dimensional tori along a resonance channel as a function of the 
rotation number (on the Poincare surface £), starting from a 
periodic orbit O a at E — 0.09. Two branches of the resonance 
channel are shown; the first branch (in blue) emerges at O a , 
the second branch (in black) appears at the bifurcation point 
of the first branch at lo/ty = 0.481422634. The points of fre- 
quency halving bifurcations "A" and "B" are the bottlenecks 
of the transition from the trapping to roaming stage. 



can be estimated as A ~ 1/AT, where N is the number 
of returns to the surface of section E before the escape. 
In Fig. [9] the proximity of the torus to a 2:5 resonance 
zone suggests the influence of a periodic orbit with or 2:5 
winding number ratio (in hollow circles). For trajectory 
close to orbit O a we have A" w 150, yielding a typical 
value of A ~ 0.007. This value is inconsistent with the 
Lyapunov exponent of the nearby resonant periodic or- 
bit (which has a Lyapunov exponent of A = 0.113 per 
return to the surface of section E), indicating that other 
structures than periodic orbits are important in describ- 
ing the capture processes. Unstable two-dimensional tori 
are indeed better candidates for the escape scenario : An 
estimate of the Lyapunov exponent in the escape stage 
is consistent with the scenario of escape along unstable 
manifolds of the resonant orbit. The local rate of tran- 
sition at the onset is estimated by the largest Lyapunov 
exponent in the family. In the case shown in Fig. [JTJ it is 
close to A = 0.06. The full picture of dynamics is compli- 
cated by existence of a family of tori with varying (and 
smaller) Lyapunov exponents. 

In Fig. [12] we represent the two families of two- 
dimensional tori considered in Fig. QT] (blue and black 
curves). First, the organizing periodic orbits O a (center 
of Fig. IT2")) and the resonance 2:5 (exterior spheres) are 
located. The projections of the two families of loops in S 
are plotted in the three dimensional space (i?i, Pi, i? 2 ). 
Meridians of the surfaces are invariant under the Poincare 
map (i.e. they are invariant loops 7). The first family 
of tori (blue curve in Fig. ITT]) starts from the central 
periodic orbit O a and continues outwards as the rota- 
tion number decreases from the value of uj^^/tt = 0.49 
(see Table [J and the blue curve in Fig. QTJ). The first 
loops of this family have zero Lyapunov exponent (the 
ones with lo/tt between 0.49 and 0.481). At the bifurca- 
tion point (uj/tv = 0.481), the second family (black curve 
on Fig. 1 1 1 [) branches off of the first one and continues nor- 
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mally (with zero Lyapunov exponent). The continuation 
of the first branch of tori is now normally hyperbolic (see 
Fig. [IT]) from uj/tt between 0.481 to 0.432, while the new 
branch of frequency halved loops is at first elliptically 
stable. The bifurcation at ui/ir = 0.481 is a frequency- 
halving, since the emerging loop winds around the orig- 
inal one twice, or in other words, has half the rotation 
number. This process is very general and we expect it 
to occur in the vicinity of any periodic orbit with sev- 
eral elliptic stability degrees of freedom. The family of 
tori has singularities at some specific rotation numbers, 
but the manifold can typically be continued across them, 
and therefore seems to be robust. The behaviour of the 
second branch of this family of tori as it approaches the 
rational rotation number lu/tt = 2/5 was investigated. A 
nontrivial foliation of invariant loops in the vicinity of 
a 2:5 periodic orbit is shown in Fig. [13J Allowed by di- 
mensional analysis, a possible scenario is that this family 
of tori is heteroclinic to the invariant manifolds of other 
invariant tori, related to periodic points in 2:5 resonance 
with O a . However, a picture of interconnected families 
of tori, permeating bulk of the entire phase space is yet 
to emerge. 

From the numerical simulations of a large assembly of 
trajectories, the following assumptions emerge : 1) the 
lowest order k = 1 resonance controls the rates of transi- 
tion from regular to chaotic dynamics, 2) the k = 1 reso- 
nance is a manifold that has a two-dimensional "back- 
bone" manifold, in analogy with resonance manifolds 
of integrable Hamiltonian systems, and 3) regular-to- 
chaotic transition occurs at the point where there is a 
transition in the normal stability of this manifold. From 
these assumptions, the typical scenario for escape after 
trapping by a weakly hyperbolic family of tori, is the fol- 
lowing one: First the trajectory evolves in a regular re- 
gion until it finds an exit channel (the transition stage) in 
the form of a manifold of normally hyperbolic invariant 
two-dimensional tori, and follows along a manifold be- 
coming more chaotic progressively, as it visits invariants 
with larger hyperbolicity (Lyapunov exponent). Eventu- 
ally it is escapes to a strongly chaotic region using the 
unstable manifold of a hyperbolic periodic orbit with a 
large Lyapunov exponent. 



IV. CONCLUSIONS 

In contrast to collinear OCS where the phase space is 
roughly divided into islands and chaotic seas, the phase 
space of planar OCS exhibits a complex ocean with cur- 
rents, reefs and shoals which slow down the progress to- 
ward energy equilibration. In this article, we have identi- 
fied these structures and their linear stability properties. 
Principal among them are two-dimensional invariant tori 
which occur in families and can be parametrized by their 
rotation numbers. These structures are organized around 
periodic orbits which provide the backbone to the dy- 
namics. By trapping trajectories temporarily, they act 
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FIG. 12: Geometry of the families of two-dimensional tori of 
Fig. [TT] (for E — 0.09). The two branches of tori are displayed 
in the (R\,P\,R,2) projection of their Poincare section (conse- 
quently one branch is composed of one-dimensional curves). 
The first branch (nearly horizontal and corresponding to the 
blue curve in Fig. Hip emerges close to the fixed point O a 
(in the center of the figure). The second branch (nearly 
vertical and corresponding to the black curve in Fig. Ill[) 
emerges at the bifurcation point of the first branch (with 
lu/tv — 0.481422634). We also represent the points (in black) 
of the 2:5 resonant periodic orbit which obstructs the con- 
tinuation of the second branch. Note that this figure links 
Fig. [9] (where two different projections are plotted and only 
one torus shown) with Fig. 1111 



as bottlenecks to the exploration of larger parts of phase 
space. Our work also makes explicit the mechanisms by 
which trajectories are trapped and by which they escape 
from the trap. 
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APPENDIX A: DISCRETE SYMMETRIES 

The time-reversibility of Hamiltonian ([1]) induces 
discrete symmetries which are taken into account to 
uniquely define invariant points on the surface of section 
£ and to to evaluate multiplicities of periodic orbits. 

Time reversibility symmetry, valid in each degree of 
freedom individually, induces "pmm" (in crystallographic 
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FIG. 13: Approach to a rational rotation number 2:5 in 
the second branch of a family of tori of periodic orbit O a , 
E = 0.09 (see also Fig. [11]) Rotation numbers, from right 
to left are 2u/n = 0.4148,0.4110,0.4104. Nontrivial foliation 
around manifolds of a resonant 2:5 periodic orbit can be seen. 
Note that this figure is related to frequency-halved loops in 
the second branch of Fig. 111! therefore we have doubled the 
original rotation numbers. 



classification) symmetry group C 2v which acts on intrin- 
sic coordinates P\, P 2 , a and P a , while R\ and R 2 are 
left invariant. Elements of Civ are identity e, reflection 
a"x, reflection a 2 , and inversion i, defined as : 

e{P 1 ,P 2 ,a,P a ) = (P u P2,a,P a ), 
ai(Pi,P 2 ,a,P Q ) = (P u P 2 ,2n-a,-P a ), 
a 2 (Px,P 2 ,a,P a ) = (-P 1 ,-P 2 ,2n-a,P a ), 
i(P u P 2 ,a,P a ) = (-P u -P 2 ,a,-P a ). 

This discrete symmetry is useful for the method of surface 
of section, because it allows to relate points x in phase 
space with P a = and P a < 0, not on the surface E, 
with points ci(x) which are on the surface E. 

In addition to exact discrete symmetries discussed 
above, the specific form of potential energy ([2]) in- 
duces an approximate Ri~-R 2 reflection symmetry as seen 
in Fig. [TJ Equation © can be written in the form of 

V = D 1 V U (R 1 ;^,R° 1 )+D 2 V M (R 2 ;I3 2> R^) 
+D 3 V M (R 3 ; fa, R° 3 ) + ViiRuRa, R 3 ), 

where V M = [1 - exp(-/3(R - R°))} 2 . Using D = (D 1 + 
D 2 )/2, and SD = (D 2 — Di)/2, the potential is rewritten 
as V(Ri,R 2 ,R 3 ) = U (Ri,R 2 ) + U I {R 1 ,R 2 ,R 3 ), where 

U a = D (V M (R i; fa, R°) + V M (R 2 ; (3 2 , R° 2 )) , 
U I = SD (V M (R 2 ;fa, R° 2 ) - V M {R 1 ; fa, i??)) 
Vi(Ri,R 2 ,R 3 ). 



-D 3 V M (R 3 ;(3 3 ,R° 3 ) 



This partition quantifies the approximate symmetry. The 
non- vanishing parameters SD, D 3 and A measure the 
deviation from the exactly symmetry. For the planar 
OCS, these parameters are SD = 0.065, D 3 = 0.16 and 
A = 0.2 compared with D = 0.348. 

With respect to linear transformations, Morse poten- 
tials transform as 

V M (aR + b; /3, R°) = V M (R; a/3, (R° - b)/a). 



Considering the linear transformations of the coordinates 
Ri and R 2 given by L(Ri,R 2 ) = (aiR 2 + b\,a 2 Ri + b 2 ), 
the symmetry line is obtained by requiring that 

U (L(R 1 ,R 2 )) = Uo(Ri,R 2 ). 

The solution is obtained in terms of parameters (3\ and 
f3 2 , and in particular f3\j(3 2 f=a 0.9, and the symmetry is 
then given by the equation : 



R 2 = ^{Rx-R\) + R%. 

P2 



In case of an exact symmetry, the symmetry line would 
be a natural boundary of the elementary cell of the dy- 
namics. All orbits could be classified with respect to this 
symmetry as having a symmetric partner, or being self 
symmetric, as usually. When the symmetry is only ap- 
proximate the cell boundary argument is no longer valid, 
but the orbits can still be classified in this way, in par- 
ticular, with regards to their degeneracy. 



APPENDIX B: METHODOLOGY: 
DETERMINATION OF INVARIANT TORI AND 
THEIR LINEAR STABILITY PROPERTIES 

We briefly summarize the method we used to compute 
two dimensional invariant tori of a Hamiltonian system. 
We have seen that this is equivalent to determining closed 
invariant curves (loops) of the Poincare map on the cho- 
sen surface of section E. Furthermore we compute the 
linear stability properties of such objects. This method 
follows the one described in Ref. 14811. 



1. Determination of invariant tori 

In order to determine two-dimensional tori, we use the 
fact that the type of internal dynamics on T 1 is likely 
to be a rotation. We assume that the Poincare map Ty, 
has an invariant curve with an irrational rotation number 
oj, and that there exists a map (at least continuous) x : 
T 1 1— > E such that Denjoy's theorem states that such 
a rotation number ui can be defined. Let (^(T 1 , E) be the 
space of continuous functions from T 1 in E, and let us 
define the linear map T w : C(T\E) i-> CfT 1 , E) as the 
translation by oj, i.e. (T w x)(#) = x(0 + uj). We define 
F : CQT^E) h+ CpP.E) as 



F(x)(0)=^ E (x(0))-(T w x)(0). 



(Bl) 



It is clear that zeros of F in Ci^f 1 , E) correspond to (con- 
tinuous) invariant curves of rotation number u>. The 
determination of two-dimensional invariant tori follows 
from the search of zeros of this functional. 

First we expand x((9) in a Fourier series with real co- 
efficients, 



x(<9) 



ao 
2 



fe>0 



(at cos irk8 + hk sin Trk8) , (B2) 
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where a k , G R™ for 6 N (n being the dimension of 
the flow) and x(#) is a periodic function with period 2, 
i.e. x(0 + 2) = x(9). We truncate these series at a fixed 
value of N, and determine an approximation to the 27V+1 
unknown coefficients ao, a*,, and for 1 < k < TV. We 
construct the discretized version of Eq. (J9j) by considering 
a mesh of 27V + 1 points on T 1 : 



0; 



2j 



27V + 1 



for < j < 27V, 



where we notice that for numerical stability reasons, 
the length of T 1 is taken as 2. Given the Fourier co- 
efficients afe, bfc, the coordinates x(0j) are expressed 
as linear functions of the coefficients a^, b^, i.e. 
x(0j) = </>({a fe },{b fc },j), given by Eq. jB2|. Accord- 
ingly, JFs(x(0j )) and Eq. ((9|) can be considered as func- 
tions of the coefficients a/j, b^: 



F,({a fc },{b fc }, 



= F s (<f>({a k },{h k },j)) 
-<j>({a k },{b k },j + i(u)), 



for < j < 27V and where = (27V + l)u)/2. The co- 
efficients a k , h k are the unknowns in the above equation. 

We solve F = using a Newton's iterative algorithm. 
At each iteration, it provides the corrections 5a k and 5h k 
to be added to the a^ and h k obtained from the previous 
iteration. We approximate (Sa, Sh) as a solution of the 
following equation: 



F 3 -(a,b,i/) 



da k 



Sa k 



dFj ru 5F, . 

— ^Sh k + — ^Sw = 0, 

oh k ouj 



where a = (ao, ai, . . . , ajv) and b = (bi, . . . , bjv). The 
iteration a' = a + Sa, b' = b + Sh and io' = ui + Sui 
converges if the initial guess is close enough to the true 
solution. The above equation requires the inversion of 
the Jacobian of Fj. From the previous definitions it is 
clear that if x(0) is a Fourier series corresponding to an 
invariant curve then, for any (f e T 1 , y(0) = x(6* + (p) 
is a different Fourier scries corresponding to the same 
invariant curve as x(0). This implies that the Jacobian 
of Fj around the invariant curve has, at least, a one- 
dimensional kernel. To solve this problem we use the 
Singular Value Decomposition. Even if Newton's algo- 
rithm has converged, we cannot claim with certainty that 



a smooth two-dimensional torus has been found. We have 
noticed that crude discretization can wash out the details 
of non-smooth curves. Sometimes doubling the number 
of points in the discretization turns a convergent case 
into a divergent one. In most cases the reliability of a so- 
lution is almost certain when testing the spectrum of the 
solution (and the norm of its eigenvectors weighted by 
the frequency, penalizing high harmonics): a smooth so- 
lution should contain a unit eigenvalue. This is why it is 
also important to monitor the linear stability properties 
of the curves we obtain numerically. 

2. Linear stability properties 

In addition to the determination of the location of the 
invariant tori, we compute their linear stability proper- 
ties to obtain information on the dynamics in its (in- 
finitesimal) neighborhood, i.e. eigenvalues and eigenvec- 
tors which give at first order an approximation to the 
invariant manifolds (stable, unstable and central) near 
the invariant curve. We consider the generalized eigen- 
value problem which amounts to finding (A,i/>) such that 



7V E (x(0))V(0) = W{0 + uo). 



(B3) 



The eigenvalues A have the following properties j48[: 1) 
A = 1 is an eigenvalue of Eq. (|B3j) : the corresponding 
eigenvector is the derivative of the loop x, 2) if A is 
an eigenvalue of Eq. (|B3[) . then A exp (2bkiTUj) is also an 
eigenvalue for any k 6 Z, 3) the closure of the set of 
eigenvalues of Eq. (|B3j) is a union of circles centered at 
the origin. 

There are two unit eigenvalues in the spectrum of 
DTyi ■ The symplectic symmetry implies that the tori are 
degenerate in the linear approximation. It implies the ex- 
istence of a family of (smooth) two-dimensional tori. As 
it is usual, we expect that this family is discontinuous and 
the discontinuities are around rational rotation numbers 
ui = m/n. Numerically, once an invariant torus with a 
specific ui is found, we simply increment the frequency 
parameter ui — > u> + 8u> and restart the search. In this 
way, we determine these families of two-dimensional tori 
parametrized by their frequency ui on the Poincare sec- 
tion. More details on the algorithm are given in Ref. [5(| • 
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